# -*- coding: utf-8 -*-
"""
Created on Thu Dec  2 10:14:38 2021

@author: zhiyinan
"""

import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import Parameters1D_casadi as par 
from Model1D_casadi_crop import *
import van_Genuchten1D_casadi as vg 
from scipy import integrate
from sklearn.preprocessing import MinMaxScaler
from casadi import *
import matplotlib.font_manager as font_manager
import time
import pandas as pd
from sklearn.linear_model import Ridge
import random

# %% Data generation - function
past = 20
future = 1

def multivariate_data(dataset, target, start_index, end_index, history_size,
                      target_size,u_shape): #, single_step=False):
    data = []
    labels = []

    start_index = start_index + history_size
    if end_index is None:
        end_index = len(dataset) - target_size
        
    for i in range(start_index, end_index):
        temp_p = dataset[i-history_size:i,:]
        for j in range(target_size-1):
        # temp_p = np.concatenate(((dataset[i-history_size:i,:]).flatten(), 
        #                            (dataset[i:i+target_size-1,:u_shape]).flatten()))
            temp = np.array([dataset[i+j,0], dataset[i,1]]).reshape((1,2))
            temp_p = np.concatenate((temp_p,temp))
        data.append(temp_p)
        labels.append((target[i:i+target_size].flatten()))
    

    return np.array(data), np.array(labels)
# %% NN layers - functions
def sigmoid(x):
    return 1/(1+np.exp(-x))

def LSTMlayer_sig(warr, uarr, barr, x_t,h_tm1,c_tm1):
    
    '''
    c_tm1 = np.array([0,0]).reshape(1,2)
    h_tm1 = np.array([0,0]).reshape(1,2)
    x_t   = np.array([1]).reshape(1,1)
    
    warr.shape = (nfeature,hunits*4)
    uarr.shape = (hunits,hunits*4)
    barr.shape = (hunits*4,)
    
    '''
    
    s_t = (mtimes(x_t,warr) + mtimes(h_tm1,uarr) + barr.reshape((1, barr.shape[0])))
    hunit = uarr.shape[0]
    i  = sigmoid(s_t[:,:hunit])
    f  = sigmoid(s_t[:,1*hunit:2*hunit])
    # _c = np.tanh(s_t[:,2*hunit:3*hunit])
    _c = sigmoid(s_t[:,2*hunit:3*hunit])
    o  = sigmoid(s_t[:,3*hunit:])
    c_t = i*_c + f*c_tm1
    # h_t = o*np.tanh(c_t)
    h_t = o*sigmoid(c_t)
    return h_t, c_t

def LSTMlayer_tanh(warr, uarr, barr, x_t,h_tm1,c_tm1):
    
    '''
    c_tm1 = np.array([0,0]).reshape(1,2)
    h_tm1 = np.array([0,0]).reshape(1,2)
    x_t   = np.array([1]).reshape(1,1)
    
    warr.shape = (nfeature,hunits*4)
    uarr.shape = (hunits,hunits*4)
    barr.shape = (hunits*4,)
    
    '''
    
    s_t = (mtimes(x_t,warr) + mtimes(h_tm1,uarr) + barr.reshape((1, barr.shape[0])))
    hunit = uarr.shape[0]
    i  = sigmoid(s_t[:,:hunit])
    f  = sigmoid(s_t[:,1*hunit:2*hunit])
    _c = np.tanh(s_t[:,2*hunit:3*hunit])
    # _c = sigmoid(s_t[:,2*hunit:3*hunit])
    o  = sigmoid(s_t[:,3*hunit:])
    c_t = i*_c + f*c_tm1
    h_t = o*np.tanh(c_t)
    # h_t = o*sigmoid(c_t)
    return h_t, c_t

# # %% Import the LSTM
# LSTM = tf.keras.models.load_model("Modeling/LSTM_5s_5s_t_10nt")
# # model_3 = model_1
# def M(x):
#     weightLSTM_1 = LSTM.layers[0].get_weights()
#     warr1, uarr1, barr1 = weightLSTM_1

#     weightLSTM_2 = LSTM.layers[1].get_weights()
#     warr2, uarr2, barr2 = weightLSTM_2
    
#     weightD = LSTM.layers[2].get_weights()
#     kernel3, bias3 = weightD
    
#     n1 = 5
#     # Initialise the model with the hidden states
#     h_t1 = SX.zeros((1, n1))
#     # Long term memory
#     C_t1 = SX.zeros((1, n1))
#     for i in range(x.shape[0]):
#         xi = x[i,:].reshape((1,2))
#         L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
#         h_t1 = vertcat(h_t1, L1[0])
#         C_t1 = vertcat(C_t1, L1[1])
    
#     n2 = 5
#     # Initialise the model with the hidden states
#     h_t2 = SX.zeros((1, n2))
#     # Long term memory
#     C_t2 = SX.zeros((1, n2))
#     for i in range(h_t1.shape[0]):
#         xi = h_t1[i,:].reshape((1,n1))
#         L2 = LSTMlayer_sig(warr2, uarr2, barr2, xi, h_t2[i,:], C_t2[i,:])
#         h_t2 = vertcat(h_t2, L2[0])
#         C_t2 = vertcat(C_t2, L2[1])
#     #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
#     # Dense layer
#     out2 = []
#     for i in range(bias3.shape[0]):
#         temp2 = mtimes(h_t2[-1,:], kernel3)[i]
#     temp2 = bias3.reshape((1, bias3.shape[0])) + mtimes(h_t2[-1,:], kernel3)
#     out2 = tanh(temp2)     
    
#     return out2

# %% Define Model 1
model_1 = tf.keras.models.load_model("Modeling/p20_f1_10s_sig_2h_012_027_5_150_e8_n")
def M1(x):
    weightLSTM_1 = model_1.layers[0].get_weights()
    warr1, uarr1, barr1 = weightLSTM_1

    weightD = model_1.layers[1].get_weights()
    kernel2, bias2 = weightD
    
    n1 = 10
    # Initialise the model with the hidden states
    h_t1 = SX.zeros((1, n1))
    # Long term memory
    C_t1 = SX.zeros((1, n1))
    for i in range(x.shape[0]):
        xi = x[i,:].reshape((1,2))
        L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
        h_t1 = vertcat(h_t1, L1[0])
        C_t1 = vertcat(C_t1, L1[1])
    
    #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
    # Dense layer
    out2 = []
    for i in range(bias2.shape[0]):
        temp2 = mtimes(h_t1[-1,:], kernel2)[i]
    temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(h_t1[-1,:], kernel2)
    out2 = sigmoid(temp2)     
    
    return out2
    
# %% Define sub model 2
model_2 = tf.keras.models.load_model('Modeling/p20_f1_10s_tanh_2h_021_032_8_90_e8_n')

def M2(x):
    weightLSTM_1 = model_2.layers[0].get_weights()
    warr1, uarr1, barr1 = weightLSTM_1

    weightD = model_2.layers[1].get_weights()
    kernel2, bias2 = weightD
    
    n1 = 10
    # Initialise the model with the hidden states
    h_t1 = SX.zeros((1, n1))
    # Long term memory
    C_t1 = SX.zeros((1, n1))
    for i in range(x.shape[0]):
        xi = x[i,:].reshape((1,2))
        L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
        h_t1 = vertcat(h_t1, L1[0])
        C_t1 = vertcat(C_t1, L1[1])
    
    #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
    # Dense layer
    out2 = []
    for i in range(bias2.shape[0]):
        temp2 = mtimes(h_t1[-1,:], kernel2)[i]
    temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(h_t1[-1,:], kernel2)
    out2 = tanh(temp2)     
    
    return out2

# %% Define sub model 3
model_3 = tf.keras.models.load_model("Modeling/p20_f1_5s_3s_sig_2h_029_038_4_35_e7_n")

def M3(x):
    weightLSTM_1 = model_3.layers[0].get_weights()
    warr1, uarr1, barr1 = weightLSTM_1

    weightLSTM_2 = model_3.layers[1].get_weights()
    warr2, uarr2, barr2 = weightLSTM_2
    
    weightD = model_3.layers[2].get_weights()
    kernel3, bias3 = weightD
    
    n1 = 5
    # Initialise the model with the hidden states
    h_t1 = SX.zeros((1, n1))
    # Long term memory
    C_t1 = SX.zeros((1, n1))
    for i in range(x.shape[0]):
        xi = x[i,:].reshape((1,2))
        L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
        h_t1 = vertcat(h_t1, L1[0])
        C_t1 = vertcat(C_t1, L1[1])
    
    n2 = 3
    # Initialise the model with the hidden states
    h_t2 = SX.zeros((1, n2))
    # Long term memory
    C_t2 = SX.zeros((1, n2))
    for i in range(h_t1.shape[0]):
        xi = h_t1[i,:].reshape((1,n1))
        L2 = LSTMlayer_sig(warr2, uarr2, barr2, xi, h_t2[i,:], C_t2[i,:])
        h_t2 = vertcat(h_t2, L2[0])
        C_t2 = vertcat(C_t2, L2[1])
    #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
    # Dense layer
    out2 = []
    for i in range(bias3.shape[0]):
        temp2 = mtimes(h_t2[-1,:], kernel3)[i]
    temp2 = bias3.reshape((1, bias3.shape[0])) + mtimes(h_t2[-1,:], kernel3)
    out2 = sigmoid(temp2)
    
    return out2

# %%
model_a = tf.keras.models.load_model("Modeling/2nd_layer_4in_5s_s_sub_noise")
def Ma(x):
    weightD1 = model_a.layers[0].get_weights()
    kernel1, bias1 = weightD1

    weightD2 = model_a.layers[1].get_weights()
    kernel2, bias2 = weightD2
    
    # weightD3 = model_a.layers[2].get_weights()
    # kernel3, bias3 = weightD3
    
    # weightD4 = model_a.layers[3].get_weights()
    # kernel4, bias4 = weightD4
    
    # n1 = 3
    # out1 = []
    # for i in range(bias1.shape[0]):
    #     temp1 = mtimes(x, kernel1)[i]
    temp1 = bias1.reshape((1, bias1.shape[0])) + mtimes(x.reshape((1,x.shape[0])), kernel1)
    out1 = sigmoid(temp1)   
    
    temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(out1, kernel2)
    out2 = sigmoid(temp2)     
    
    # temp3 = bias3.reshape((1, bias3.shape[0])) + mtimes(out2, kernel3)
    # out3 = tanh(temp3)  
    
    # temp4 = bias4.reshape((1, bias4.shape[0])) + mtimes(out3, kernel4)
    # out4 = sigmoid(temp4)  
    
    return out2


# %% Richards Eqn Simulation with disturbances

past = 20
future = 1

y_index = 19

soilPars = loamySoil()

h0 = -0.96
x_i = np.ones(26)*h0

lbu = -1e-06  #  unit [m/s]
ubu = 0  #  unit [m/s]
lby = vg.volumetricMoisture(-80)
uby = vg.volumetricMoisture(-0.1)

NN_past = np.loadtxt("NN_in_0.txt")
NN_past[0,0] = 0
y0 = NN_past[0,1]
NN_past[:,0] = (NN_past[:,0] - lbu)/(ubu-lbu)
NN_past[:,1] = (NN_past[:,1] - lby)/(uby-lby)
# Define the original model for optimal 
Ny = 1
totalDepth, axialNodes, axialDistance, totalNodes = par.spatialVariables_1D()  #Spatial parameters
samplingTime, samplingTimeInternal, internalTimeSteps, totTimeSteps = \
                                 par.temporalVariable_dataGen(1200*7200)

def ODE(head, irrigAmount, cropCoeff, refEvap):
    origin = Richards1D_casadi(head, irrigAmount, cropCoeff, refEvap)
    error = DM.zeros(26)
    # for i in range(26):
    #     error[i] = np.random.choice([-1,1]) * 0.1*np.random.random_sample()*origin[i]
    return origin + error

t = MX.sym('t')
x = MX.sym('x', totalNodes)
I = MX.sym('I')
C = MX.sym('C')
rE = MX.sym('rE')

p = vertcat(I, C, rE)
ode = {'x': x, 'p': p, 'ode':ODE(x, p[0], p[1], p[2])}
opts = {'tf': samplingTimeInternal, 'regularity_check': True}
F = integrator('F', 'cvodes', ode, opts)

def sim_2h_MX(head, irrigAmount, cropCoeff, refEvap):
    for t in range(internalTimeSteps):
        r = F(x0 = head, p = vertcat(irrigAmount, cropCoeff, refEvap))
        head = r["xf"]
        # y = vg.volumetricMoisture(x_i[y_index])
    return head

cropCoeff = 0.88
refEvap = 3.102e-08


#%%  Import the weather condition 
# data = pd.read_csv('C:/Users/zhiyinan/OneDrive - ualberta.ca/Irrigation/crop_model_1D\Weather.csv')
# refEvap = (data.iloc[10:20,5].values)/1e3/(24*60*60)
random.seed(5)
np.random.seed(5)

# Nd = int(totTimeSteps/12)
# refEvap = np.random.normal(3, 0.2, Nd)/1e3/(24*60*60)
# Evap = np.zeros(totTimeSteps)
# for i in range(Nd):
#     Evap[i*12:i*12+12] = refEvap[i]
# np.savetxt("C:/Users/zhiyinan/OneDrive - ualberta.ca/Irrigation/crop_model_1D/ZMPC Simulation Results/Evap.txt", Evap)
# data_hr = pd.read_csv('C:/Users/zhiyinan/OneDrive - ualberta.ca/Irrigation/crop_model_1D\Rain_hourly.csv')
# precip = data_hr.iloc[240:480,2].values/1e3/(60*60*24)


Precip = np.zeros((Nd,12))

# for i in range(Nd):
#     start = random.randrange(12)
#     end = random.randrange(start, 12)
#     for j in range(start, end):
#         Precip[i, j] = -(random.random()*(20-5)+5)/1e3/(24*60*60)

Precip = Precip.flatten()
# Precip = np.zeros(int(Nsim/Delta))

# for i in range(int(Nsim/Delta)):
#     Precip[i] = random.choice([0,0,0,0,1]) * (random.random()*(2-0.2)+0.2)/1e3/(24*60*60)

# plt.figure(8)
# plt.step(np.arange(100), Precip[:100]) # np.arange(Precip.shape[0])
# np.savetxt("C:/Users/zhiyinan/OneDrive - ualberta.ca/Irrigation/crop_model_1D/ZMPC Simulation Results/Precip.txt", Precip)
# lzone_store = DM.zeros((int(Nsim/Delta),Ny))
# uzone_store = DM.zeros((int(Nsim/Delta),Ny))
# start = time.time()

u_sim = np.zeros((totTimeSteps,1))

i=0
while i < totTimeSteps:
    irr_freq = int(50*60*60/samplingTime)
    ulist = np.array([-7e-6, -6e-6, -5e-6, -4e-6, -3e-6, -2e-6, -1.5e-6, -1e-6, -0.5e-6])
    # ulist = [-2/1000/3600, -1.5/1000/3600, -0.9/1000/3600, -1.2/1000/3600, -0.5/1000/3600]
    t = irr_freq*random.randrange(1,3) + random.randrange(irr_freq)
    u_sim[i] = random.choice(ulist)
    i = i+t
    print(i)

x_sim = np.zeros((totTimeSteps+1, totalNodes))
x_sim[0,:] = x_i
y_sim = np.zeros((totTimeSteps+1, Ny))
y_sim[0,:] = vg.volumetricMoisture(x_i[y_index])
for t in range(totTimeSteps):
    x_sim[t+1,:] = np.array(sim_2h_MX(x_sim[t,:], u_sim[t,:], cropCoeff, refEvap)).flatten()
    y_sim[t+1,:] = vg.volumetricMoisture(x_sim[t+1, y_index])
# %% Load the data
y_index = 19
Ny = 1
uv = np.loadtxt("Modeling/u_Richards_bias.txt")
yv = np.loadtxt("Modeling/y_Richards_bias.txt")

t_pred = 20

np.random.seed(5)
yn = 0.1*np.random.choice([-1,1],(yv.shape))*np.random.random((yv.shape))
yv_n = yv*(1+yn)

y_test = yv[1:]
x_test = np.zeros((uv.shape[0],2))
x_test[:,0] = uv
x_test[:,1] = yv_n[:-1]


scaler_xt = MinMaxScaler()
x_test = scaler_xt.fit_transform(x_test)
y_test = y_test.reshape((y_test.shape[0],1))
scaler_yt = MinMaxScaler()
y_test = scaler_yt.fit_transform(y_test)

x_init, y_t = multivariate_data(x_test, y_test, 0, None, past, future,1)
# yo = scaler_yt.inverse_transform(y_t.reshape((y_t.shape[0], y_t.shape[1])))

OUTPUT_SIZE = future*Ny

# x_init = []
y_act = []
# for i in range(x_t.shape[0]):
#     if (i+1)%t_pred ==1:
#         x_init.append(x_t[i,:,:])
#         y_act.append(y_t[i:i+t_pred,:].flatten())
for i in range(x_t.shape[0]-t_pred):
    y_act.append(y_t[i:i+t_pred,:].flatten())
# x_init = np.array(x_init[:-1])
y_act = np.array(y_act[:-1])


# %% Test the bias policy - constant

error = []
y_save = np.zeros((x_init.shape[0], t_pred))
b_save = []

e = 0
b = 0
update = 1
n=500
for t in range(n): 
    print(t)
    y_sim_1 = np.zeros((t_pred,1))
    y_NN = np.zeros((t_pred,1))
    NN0_1 = x_init[t]
    y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M3(NN0_1)), NN0_1[-1, 1]])  #, DM(M4(NN0_1)), DM(M5(NN0_1))])
    y_sim_1[0,:] = scaler_yt.inverse_transform(np.array(DM(Ma(y_sub)))).flatten() + b
    y_NN[0] = scaler_yt.transform(y_sim_1[0].reshape((1,1))).flatten() 
    e += y_act[t, 0] - y_sim_1[0,0]
    u_sim = x_t[t:t+t_pred,-1,0]
    for i in range(t_pred-1):
        NN0_1 = np.concatenate((NN0_1[1:], np.array([u_sim[i], y_NN[i,0]]).reshape((1,2))))
        y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M3(NN0_1)), NN0_1[-1, 1]])  #, DM(M4(NN0_1)), DM(M5(NN0_1))])
        y_sim_1[i+1,:] = scaler_yt.inverse_transform(np.array(DM(Ma(y_sub)))).flatten() + b
        y_NN[i+1] = scaler_yt.transform(y_sim_1[i+1].reshape((1,1))).flatten()
        e += y_act[t, i+1] - y_sim_1[i+1,0]
    # error.append(e)
    # e = 0
    
        if (i+2)%update == 0:
            # b = e/update  
            b = sign(e)*min(abs(e/update), 0.2)
            b_save.append(b)
            error.append(e)
            e = 0
    # yp_1 = scaler_yt.inverse_transform(y_sim_1) + b
    y_save[t,:] = y_sim_1.flatten()

diff = y_act[:n,:] - y_save[:n,:]
print(sum(sum(abs(diff)))/n/y_act.shape[1])
plt.figure(1)
plt.plot(b_save) # 
plt.plot(np.array(error)/update)

plt.figure(2)
plt.plot(diff[:,0])
# %%

def error_model(y_sim, y_pred):
    W = []
    W0 = []
    W_lb = []
    W_ub = []
    f = 0
    # g = []
    # g_lb = []
    # g_ub = []
    
    coef = SX.sym('coef', 2)
    W = vertcat(W, coef[0])
    W0 = vertcat(W0, 1)
    W_lb = vertcat(W_lb, 0.8)
    W_ub = vertcat(W_ub, 1.5)
    
    W = vertcat(W, coef[1])
    W0 = vertcat(W0, 0)
    W_lb = vertcat(W_lb, -0.2)
    W_ub = vertcat(W_ub, 0.3)
    n = y_pred.shape[0]
    for i in range(n):
        f += (y_sim[i] - (y_pred[i]*coef[0]+coef[1]))**2
    nlp = {'x':W,'f':f}
    S = nlpsol('S','ipopt',nlp)
    r = S(x0 = W0, lbx = W_lb, ubx = W_ub) 
    return r['x'][0], r['x'][1]


error = []
y_save = np.zeros((x_init.shape[0], t_pred))
a_save = []
b_save = []

e = 0
# a = 1
# b = 0
update = 10
for t in range(n): 
    print(t)
    a = 1
    b = 0
    y_sim_1 = np.zeros((t_pred,1))
    y_NN = np.zeros((t_pred,1))
    NN0_1 = x_init[t]
    y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M3(NN0_1)), NN0_1[-1, 1]])  #, DM(M4(NN0_1)), DM(M5(NN0_1))])
    y_sim_1[0,:] = scaler_yt.inverse_transform(np.array(DM(Ma(y_sub)))).flatten()*a + b
    y_NN[0] = scaler_yt.transform(y_sim_1[0].reshape((1,1))).flatten() 
    e += y_act[t, 0] - y_sim_1[0,0]
    u_sim = x_t[t:t+t_pred,-1,0]
    for i in range(t_pred-1):
        NN0_1 = np.concatenate((NN0_1[1:], np.array([u_sim[i], y_NN[i,0]]).reshape((1,2))))
        y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M3(NN0_1)), NN0_1[-1, 1]])  #, DM(M4(NN0_1)), DM(M5(NN0_1))])
        y_sim_1[i+1,:] = scaler_yt.inverse_transform(np.array(DM(Ma(y_sub)))).flatten()*a + b
        y_NN[i+1] = scaler_yt.transform(y_sim_1[i+1].reshape((1,1))).flatten()
        e += y_act[t, i+1] - y_sim_1[i+1, 0]
        if (i+2)%update == 0:
            a,b = error_model(y_act[t, i-update+1:i+1], y_sim_1[i-update+1:i+1,0])
            a_save.append(a)
            b_save.append(b)
            error.append(e)
            e = 0
    # yp_1 = scaler_yt.inverse_transform(y_sim_1) + b
    y_save[t,:] = y_sim_1.flatten()

diff = y_act[:n,:] - y_save[:n,:]
print(sum(sum(abs(diff)))/n/y_act.shape[1])

# plt.figure(1)
# plt.plot(np.array(error)/update)

# %% Bias Design
a = 900
uf = 10 # Update frequency of the bias
yp_LSTM = np.zeros((a, 1))
yp_2l = np.zeros((a, 1))
for t in range(a):
    print(t) 
    y_sim_LSTM = np.zeros((t_pred,1))
    y_sim_2l = np.zeros((t_pred,1))
    
    NN0_1 = x_t[t,:]
    u_sim = x_t[t:t+t_pred,-1,0]
    
    y_sim_LSTM[0,:] = np.array(DM(M(NN0_1))).flatten()
    y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M3(NN0_1)), NN0_1[-1,1]])  #, DM(M4(NN0_1)), DM(M5(NN0_1))])
    y_sim_2l[0,:] = np.array(Ma(y_sub)).flatten()
    for i in range(t_pred-1):
        NN0_1 = np.concatenate((NN0_1[1:], np.array([u_sim[i], y_sim_1[i,0]]).reshape((1,2))))
        y_sim_LSTM[i+1,:] = np.array(DM(M(NN0_1))).flatten()
        
        y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M3(NN0_1)), NN0_1[-1,1]])  
        y_sim_2l[i+1,:] = np.array(Ma(y_sub)).flatten()
    
    y_LSTM = scaler_yt.inverse_transform(y_sim_LSTM).flatten()
    y_2l = scaler_yt.inverse_transform(y_sim_2l).flatten()
    yp_LSTM[t,:] = y_LSTM
    error[t] = sum((y_act[t:t+t_pred].flatten() - yp_1)**2/y_act[t:t+t_pred].flatten())
    
    

    for i in range(t_pred-1):
        NN0_1 = np.concatenate((NN0_1[1:], np.array([u_sim[i], y_sim_1[i,0]]).reshape((1,2))))
        y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M3(NN0_1)), NN0_1[-1,1]])  
        y_sim_1[i+1,:] = np.array(Ma(y_sub)).flatten()
    yp_1 = scaler_yt.inverse_transform(y_sim_1)
    y_save[t0,:] = yp_1.flatten()
    error[t0] = sum((y_act[t0:t0+t_pred] - yp_1)**2/y_act[t0:t0+t_pred])
    
    















